Aplicaciones del Big Data en la Empresa 02
Ivonne Yañez Mendoza
Master Big Data & Business Analytics, Universidad Complutense de MadridOctubre 2022
Preparación del entorno de trabajo y carga de librerÃas
# Library loading
rm(list = ls())
suppressPackageStartupMessages({
library(data.table)
library(dplyr)
library(caret)
library(scales)
library(ggplot2)
library(stringi)
library(stringr)
library(dataPreparation)
library(knitr)
library(kableExtra)
library(ggpubr)
library(tictoc)
library(ggeasy)
library(lubridate)
library(inspectdf)
library(ranger)
library(gbm)
})1 Introducción al problema
Utilizando los datos provistos por el Ministerio del agua de Tanzania, se requiere construir un modelo que sea capaz de predecir cuales bombas de agua están operativas, operativas pero que necesitan reparación o están dañadas, basadas en un set de datos train.
En el primer modelo solo se realizaron pruebas con las variables numéricas, obteniendo un scoring aceptable pero mejorable.
En este segundo modelo se realizaran combinaciones de variables y transformaciones buscando mejoras y optimizaciones.
2 Exploración inicial al conjunto de datos
De forma original se incluye los dataset train y test mas el archivo con las labels.
# Data Loading
#Fichero datos train
vtrain <- fread("train_set.csv")
vtrain$flag <- 1 # Columna que indica si es parte del set train (1) test (0)
#Fichero datos test
vtest <- fread("test_set.csv")
vtest$flag <- 0 # Columna que indica si es parte del set train (1) test (0)
#Fichero con labels o objetivo
vlabels <-fread("labels.csv")
# Se unen las labels con el set de datos train
train <- merge(vlabels, vtrain)
# Se unen ambos datasets (train y test) lo cual es lo recomendado para mas adelante trabajar en Feature engineering
datos <- as.data.table(rbind(vtrain, vtest))
#Comprobacion
head(datos)head(vlabels)#Distribucion de los datos
str(datos)## Classes 'data.table' and 'data.frame': 74250 obs. of 41 variables:
## $ id : int 69572 8776 34310 67743 19728 9944 19816 54551 53934 46144 ...
## $ amount_tsh : num 6000 0 25 0 0 20 0 0 0 0 ...
## $ date_recorded : IDate, format: "2011-03-14" "2013-03-06" ...
## $ funder : chr "Roman" "Grumeti" "Lottery Club" "Unicef" ...
## $ gps_height : int 1390 1399 686 263 0 0 0 0 0 0 ...
## $ installer : chr "Roman" "GRUMETI" "World vision" "UNICEF" ...
## $ longitude : num 34.9 34.7 37.5 38.5 31.1 ...
## $ latitude : num -9.86 -2.15 -3.82 -11.16 -1.83 ...
## $ wpt_name : chr "none" "Zahanati" "Kwa Mahundi" "Zahanati Ya Nanyumbu" ...
## $ num_private : int 0 0 0 0 0 0 0 0 0 0 ...
## $ basin : chr "Lake Nyasa" "Lake Victoria" "Pangani" "Ruvuma / Southern Coast" ...
## $ subvillage : chr "Mnyusi B" "Nyamara" "Majengo" "Mahakamani" ...
## $ region : chr "Iringa" "Mara" "Manyara" "Mtwara" ...
## $ region_code : int 11 20 21 90 18 4 17 17 14 18 ...
## $ district_code : int 5 2 4 63 1 8 3 3 6 1 ...
## $ lga : chr "Ludewa" "Serengeti" "Simanjiro" "Nanyumbu" ...
## $ ward : chr "Mundindi" "Natta" "Ngorika" "Nanyumbu" ...
## $ population : int 109 280 250 58 0 1 0 0 0 0 ...
## $ public_meeting : logi TRUE NA TRUE TRUE TRUE TRUE ...
## $ recorded_by : chr "GeoData Consultants Ltd" "GeoData Consultants Ltd" "GeoData Consultants Ltd" "GeoData Consultants Ltd" ...
## $ scheme_management : chr "VWC" "Other" "VWC" "VWC" ...
## $ scheme_name : chr "Roman" "" "Nyumba ya mungu pipe scheme" "" ...
## $ permit : logi FALSE TRUE TRUE TRUE TRUE TRUE ...
## $ construction_year : int 1999 2010 2009 1986 0 2009 0 0 0 0 ...
## $ extraction_type : chr "gravity" "gravity" "gravity" "submersible" ...
## $ extraction_type_group: chr "gravity" "gravity" "gravity" "submersible" ...
## $ extraction_type_class: chr "gravity" "gravity" "gravity" "submersible" ...
## $ management : chr "vwc" "wug" "vwc" "vwc" ...
## $ management_group : chr "user-group" "user-group" "user-group" "user-group" ...
## $ payment : chr "pay annually" "never pay" "pay per bucket" "never pay" ...
## $ payment_type : chr "annually" "never pay" "per bucket" "never pay" ...
## $ water_quality : chr "soft" "soft" "soft" "soft" ...
## $ quality_group : chr "good" "good" "good" "good" ...
## $ quantity : chr "enough" "insufficient" "enough" "dry" ...
## $ quantity_group : chr "enough" "insufficient" "enough" "dry" ...
## $ source : chr "spring" "rainwater harvesting" "dam" "machine dbh" ...
## $ source_type : chr "spring" "rainwater harvesting" "dam" "borehole" ...
## $ source_class : chr "groundwater" "surface" "surface" "groundwater" ...
## $ waterpoint_type : chr "communal standpipe" "communal standpipe" "communal standpipe multiple" "communal standpipe multiple" ...
## $ waterpoint_type_group: chr "communal standpipe" "communal standpipe" "communal standpipe" "communal standpipe" ...
## $ flag : num 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, ".internal.selfref")=<externalptr>Se despliega la distribución de la variable objetivo (contenida en labels) Se observa algo de desbalance el cual no deberÃa afectar demasiado.
cuenta <- vlabels %>% count(status_group)
porcentaje <- round( prop.table(table(vlabels$status_group))*100, 2)
kable(cuenta, col.names = c('status_group', 'count'))| status_group | count |
|---|---|
| functional | 32259 |
| functional needs repair | 4317 |
| non functional | 22824 |
kable(porcentaje, col.names = c('status_group', '%'))| status_group | % |
|---|---|
| functional | 54.31 |
| functional needs repair | 7.27 |
| non functional | 38.42 |
barplot(porcentaje, col=rgb(0.2,0.4,0.6,0.6))
3 Exploratory data analysis EDA
# categorical plot
x <- inspect_cat(datos)
show_plot(x)
# correlations in numeric columns
x <- inspect_cor(datos)## Warning: Columns with 0 variance found: flagshow_plot(x)
# feature imbalance bar plot
x <- inspect_imb(datos)
show_plot(x)
# memory usage barplot
x <- inspect_mem(datos)
show_plot(x)
# missingness barplot
x <- inspect_na(datos)
show_plot(x)
# histograms for numeric columns
x <- inspect_num(datos)
show_plot(x)
# barplot of column types
x <- inspect_types(datos)
show_plot(x)
3.1 Observaciones EDA
Una gran parte de las features son de tipo categóricas.
Se deben explorar los missing values (en gris) en el primer gráfico de frecuencias de las categóricas.
wpt_name, subvillage, scheme_name e installer son las que ocupan mas espacio en memoria, si se revisan en el gráfico de frecuencias de las categóricas, son las que mas categorÃas contienen por cada una de las variables.
public_meeting y permit muestran columnas con % de NA (5.6% y 5.1% respectivamente)
La feature construction_year tiene valores en 0, es decir, no hay información codificada del año de construcción de la bomba de agua.
También se deben explorar las demás variables numéricas con 0 para determinar que hacer con ellas.
4 Feature engineering
En el modelo anterior se trabajaron solo con variables numéricas, en este modelo se exploraran numéricas y categóricas.
Se exploran algunas columnas como categóricas para construir el modelo:
# Remueve atributos que no se usaran
datos$num_private <- NULL
datos$amount_tsh <- NULL
datos$wpt_name <- NULL
datos$subvillage <- NULL
datos$region_code <- NULL
datos$lga <- NULL
datos$ward <- NULL
datos$recorded_by <- NULL
datos$scheme_name <- NULL
datos$permit <- NULL
datos$extraction_type <- NULL
datos$extraction_type_class <- NULL
datos$management_group <- NULL
datos$quality_group <- NULL
datos$quantity_group <- NULL
datos$waterpoint_type_group <- NULL
datos$source_type <- NULL
datos$date_recorded <- NULL
datos$district_code <- NULL
datos$payment_type <- NULLSe inspecciona el nuevo dataset
x <- inspect_types(datos)
show_plot(x)
Revisión de valores NA o ceros
x <- inspect_na(datos)
show_plot(x)
En este caso public_meeting es una variable de tipo booleano que contiene un 5.6% de valores en NA.
table(datos$public_meeting)##
## FALSE TRUE
## 6346 63749Se imputan los NA a la opción mayoritaria que es TRUE
datos$fe_public_meeting<- ifelse(is.na(datos$public_meeting), FALSE, datos$public_meeting)
datos$fe_public_meeting_missing<-ifelse(is.na(datos$public_meeting), 1,0)
# Eliminacion de variable original
datos$public_meeting <- NULL# Syntax
temp <- datos %>% mutate_at(vars(-c(flag)), ~na_if(., 0))
#
x <- inspect_na(temp)
show_plot(x)
Se convierten los na anteriores en ceros otra vez (para no tener problemas al construir el modelo con ranger.
temp[is.na(temp)] <- 0
head(temp)Se construyen algunas variables nuevas que pueden ser útiles a la hora de modelar. Se agregan con una fe_ para identificarlas.
La variable construction_year tiene un % importante de valores missing. Se crea una columna nueva con valores missing imputados segun la media.
temp$fe_construction_year<-round(ifelse(temp$construction_year==0, mean(temp$construction_year[temp$construction_year>0]),temp$construction_year), 0)#Combinacion de latitud y longitud
temp$fe_lonlat <- sqrt(temp$longitude^2 + temp$latitude^2)Otra variable que calcula la antigüedad de la bomba basada en su fecha de construcción
year <- year(now())
temp$fe_antiguedad <- (year - temp$fe_construction_year)Se elimina la variable construction_year
temp$construction_year <- NULLSe estudian las categóricas, hay features con un numero alto de categorÃas
# Categoricas
categoricas <- names(temp[, which(sapply(temp, is.character)), with = FALSE])
#-Frecuencias
freq_inicial <- apply(temp[, ..categoricas], 2, function(x) length(unique(x)))
freq_inicial## funder installer basin
## 2141 2411 9
## region scheme_management extraction_type_group
## 21 13 13
## management payment water_quality
## 12 7 8
## quantity source source_class
## 5 10 3
## waterpoint_type
## 7Para estos casos se decide utilizar la técnica de sustituir cada categorÃa según su frecuencia de aparición
feature: funder
cfunder<-unique(temp[ , .(.N), by = .(funder)][order(-N)])
cfunder$perc<-cfunder$N/length(temp$funder)*100
temp[ , fe_funder := .N , by = .(funder)]
temp$funder <- NULL # elimina featurefeature: installer
cinstaller<-unique(temp[ , .(.N), by = .(installer)][order(-N)])
cinstaller$perc<-cinstaller$N/length(temp$installer)*100
temp[ ,fe_installer := .N , by = .(installer)]
temp$installer <- NULL # elimina featurefeature: basin
cbasin<-unique(temp[ , .(.N), by = .(basin)][order(-N)])
cbasin$perc<-cbasin$N/length(temp$basin)*100
temp[ ,fe_basin := .N , by = .(basin)]
temp$basin <- NULL # elimina featurefeature: region
cregion<-unique(temp[ , .(.N), by = .(region)][order(-N)])
cregion$perc<-cregion$N/length(temp$region)*100
temp[ ,fe_region := .N , by = .(region)]
temp$region <- NULL # elimina featurefeature: scheme_management
cscheme<-unique(temp[ , .(.N), by = .(scheme_management)][order(-N)])
cscheme$perc<-cscheme$N/length(temp$scheme_management)*100
temp[ ,fe_scheme := .N , by = .(scheme_management)]
temp$scheme_management <- NULL # elimina featurefeature: extraction_type_group
cextraction<-unique(temp[ , .(.N), by = .(extraction_type_group)][order(-N)])
cextraction$perc<-cextraction$N/length(temp$extraction_type_group)*100
temp[ ,fe_extract := .N , by = .(extraction_type_group)]
temp$extraction_type_group <- NULL # elimina featurefeature: management
cmanagement<-unique(temp[ , .(.N), by = .(management)][order(-N)])
cmanagement$perc<-cmanagement$N/length(temp$management)*100
temp[ ,fe_management := .N , by = .(management)]
temp$management <- NULL # elimina featurefeature: payment
cpayment<-unique(temp[ , .(.N), by = .(payment)][order(-N)])
cpayment$perc<-cpayment$N/length(temp$payment)*100
temp[ ,fe_payment := .N , by = .(payment)]
temp$payment <- NULL # elimina featurefeature: water_quality
cwater<-unique(temp[ , .(.N), by = .(water_quality)][order(-N)])
cwater$perc<-cwater$N/length(temp$cwater)*100
temp[ ,fe_water_quality := .N , by = .(water_quality)]
temp$water_quality <- NULL # elimina featurefeature: quantity
cquantity<-unique(temp[ , .(.N), by = .(quantity)][order(-N)])
cquantity$perc<-cquantity$N/length(temp$quantity)*100
temp[ ,fe_quantity := .N , by = .(quantity)]
temp$quantity <- NULL # elimina featurefeature: source
csource<-unique(temp[ , .(.N), by = .(source)][order(-N)])
csource$perc<-csource$N/length(temp$source)*100
temp[ ,fe_source := .N , by = .(source)]
temp$source <- NULL # elimina featurefeature: source_class
csource_c<-unique(temp[ , .(.N), by = .(source_class)][order(-N)])
csource_c$perc<-csource_c$N/length(temp$source_class)*100
temp[ ,fe_source_class := .N , by = .(source_class)]
temp$source_class <- NULL # elimina featurefeature: waterpoint_type
cwaterpoint<-unique(temp[ , .(.N), by = .(waterpoint_type)][order(-N)])
cwaterpoint$perc<-cwaterpoint$N/length(temp$waterpoint_type)*100
temp[ ,fe_waterpoint_type := .N , by = .(waterpoint_type)]
temp$waterpoint_type <- NULL # elimina feature5 Pasos previos a modelizar
# Separa train y test segun su flag
trainset<-temp[temp$flag==1,]
#table(trainset$flag) comprobacion
testset<-temp[temp$flag==0,]
# Se combina el train set con la variable objetivo
trainset <- merge(trainset, vlabels, by ='id', sort = FALSE)
# Elimina la columna flag y la columna id
trainset$flag <- NULL
testset$flag <- NULL
trainset$id <-NULLLa columna status_group indica en palabras si es funcional o no. Se recodifica para simplificar
trainset = trainset %>%
mutate(status_group = ifelse(status_group== "functional", 0,
ifelse(status_group == "functional needs repair",1 , 2)))
table(trainset$status_group)##
## 0 1 2
## 32259 4317 228246 Construccion de modelo
6.1 Random forest con ranger
Tomando como base lo anterior y considerando las variables seleccionadas anteriormente se construye el segundo modelo.
fit <- ranger(status_group ~. ,
data = trainset,
num.trees = 300,
mtry=6,
importance = 'impurity',
write.forest = TRUE,
min.node.size = 1,
splitrule = "gini",
verbose = TRUE,
classification = TRUE,
seed=1234
)Se despliegan los resultados
print(fit)## Ranger result
##
## Call:
## ranger(status_group ~ ., data = trainset, num.trees = 300, mtry = 6, importance = "impurity", write.forest = TRUE, min.node.size = 1, splitrule = "gini", verbose = TRUE, classification = TRUE, seed = 1234)
##
## Type: Classification
## Number of trees: 300
## Sample size: 59400
## Number of independent variables: 22
## Mtry: 6
## Target node size: 1
## Variable importance mode: impurity
## Splitrule: gini
## OOB prediction error: 18.60 %Al contrario del primer modelo solo con variables numéricas, aquà el modelo2 presenta una mejorÃa considerable. Se mira la matriz de confusión para trainset.
predictions_train <- predict(fit, data = trainset)
confusionMatrix(table( trainset$status_group, predictions_train$predictions))## Confusion Matrix and Statistics
##
##
## 0 1 2
## 0 32082 85 92
## 1 390 3869 58
## 2 421 53 22350
##
## Overall Statistics
##
## Accuracy : 0.9815
## 95% CI : (0.9804, 0.9826)
## No Information Rate : 0.5538
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.9663
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Statistics by Class:
##
## Class: 0 Class: 1 Class: 2
## Sensitivity 0.9753 0.96556 0.9933
## Specificity 0.9933 0.99191 0.9872
## Pos Pred Value 0.9945 0.89622 0.9792
## Neg Pred Value 0.9701 0.99749 0.9959
## Prevalence 0.5538 0.06746 0.3788
## Detection Rate 0.5401 0.06513 0.3763
## Detection Prevalence 0.5431 0.07268 0.3842
## Balanced Accuracy 0.9843 0.97874 0.9902Se predice sobre los datos del concurso
predictions_concurso <- predict(fit, data = testset)
resultados_concurso<- as.data.frame(cbind( testset$id,predictions_concurso$predictions))
names(resultados_concurso)<-c("id", "status_group")
resultados_concurso$status_group<-ifelse(resultados_concurso$status_group==0,"functional",
ifelse(resultados_concurso$status_group==1,"functional needs repair","non functional"))
#Se guarda en el fichero que se subirá a la plataforma -->
fwrite(resultados_concurso, file = "results_model2.csv")Variables importantes
vars_imp <- fit$variable.importance
vars_imp <- as.data.frame(vars_imp)
vars_imp$myvar <- rownames(vars_imp)
vars_imp <- as.data.table(vars_imp)
setorder(vars_imp, -vars_imp)Plot de variables mas importantes
ggbarplot(vars_imp,
x = "myvar", y = "vars_imp",
#fill = 'myvar',
color = "blue",
palette = "jco",
sort.val = "asc",
sort.by.groups = FALSE,
x.text.angle = 90,
ylab = "Importancia",
xlab = 'Variable',
#legend.title = "MPG Group",
rotate = TRUE,
ggtheme = theme_minimal()
)
7 Conclusiones
Las features con mas peso predictivo para este modelo son la recategorizacion fe_quantity (cuanta cantidad de agua tiene la bomba) seguida de latitud, fe_lonlat, longitud (mismas que aparecen liderando el peso predictivo en el modelo1) ademas de fe_waterpoint_type. Es decir las recategorizaciones de categoricas segun frecuencias, se vieron reflejadas en este modelo.
Este segundo modelo tiene un scoring de 0.8172 lo cual es una mejora significativa para esta nueva predicción con respecto al modelo básico anterior. Se creará un tercer modelo buscando si se puede obtener alguna mejora adicional.
